#include "App.h"
#include "PCGSolverTest.h"
#include <sstream>

App::App()
{
}

App::~App()
{
	for (unsigned i=0; i<2;i++)
		delete _levelSetTexture[i];
	//delete _levelSetFramebuffer;
	delete _levelSetProgram;
	delete _levelSetTransferTexture;
	delete _zeroOneVertexArray;
	delete _zeroOnePositionBuffer;
	delete _gridBuffer;
	delete _gridVertexArray;
	delete _pointGridBuffer;
	delete _pointGridVertexArray;
	for (unsigned i=0; i<2; i++)
		delete _velocityTexture[i];
	for (unsigned i=0;i<2;i++)
		delete _jumpFloodingTexture[i];
	delete _gridProgram;
	delete _velocityProgram;
	delete _markFluidNodesProgram;
	delete _initJumpFloodingProgram;
	delete _jumpFloodingProgram;
	delete _jumpFloodingFinalProgram;
	delete _debugTextureProgram;
	delete _advectAndForceProgram;
	delete _levelSetMarkBoundaryProgram;

	delete _pcgSolver;
	delete _fontRenderer;
}

glm::vec2 App::cursorToTex(glm::vec2 cursor)
{
	return glm::vec2(cursor.x / _windowWidth, 1.f - cursor.y / _windowHeight);
}

void App::buttonPressed(int button, int mods)
{
	if (button == 0)
	{
		_cameraCenterStart = _cameraCenter;
		double xpos, ypos;
		glfwGetCursorPos(_window, &xpos, &ypos);
		_mousePositionStart = cursorToTex(glm::vec2(xpos, ypos));
		_cameraMoving = true;
	}
}

void App::buttonReleased(int button, int mods)
{
	if (button == 0)
	{
		_cameraMoving = false;
	}
}

void App::keyPressed(int key, int scancode, int mods)
{
	if (key == GLFW_KEY_ESCAPE)
		glfwSetWindowShouldClose(_window, GL_TRUE);
	if (key == GLFW_KEY_G)
		_drawGrid = !_drawGrid;
	if (key == GLFW_KEY_V)
		_drawVelocities = !_drawVelocities;
	if (key == GLFW_KEY_F)
		_drawFluidNodes = !_drawFluidNodes;
	if (key == GLFW_KEY_D)
		_drawDebugTexture = !_drawDebugTexture;
	if (key == GLFW_KEY_SPACE)
		_simulationRunning = !_simulationRunning;
	if (key == GLFW_KEY_RIGHT)
		_oneStepForward = true;
	if (key == GLFW_KEY_R)
		_reinitLevelSetNow = true;
}

void App::keyReleased(int key, int scancode, int mods)
{
}

void App::wheelChanged(double x, double y)
{
	if (!_cameraMoving)
	{
		_zoomFactor = _zoomFactor*powf(0.9f, (float)-y);
		updateProjMatrix();
	}
}

void App::framebufferResized(int width, int height)
{
	_windowWidth = width;
	_windowHeight = height;
	_windowRatio = width / (float)height;
	_windowResized = true;
}

int App::init(GLFWwindow* window, InputHandler* inputHandler)
{
	_window = window;
	glfwGetFramebufferSize(window, &_windowWidth, &_windowHeight);
	_windowRatio = _windowWidth / (float)_windowHeight;
	_windowResized = false;
	_zoomFactor = 1.f;
	_cameraMoving = false;
	_drawGrid = false;
	_drawVelocities = false;
	_drawFluidNodes = false;
	_drawDebugTexture = false;
	_jumpFloodingSrc = 0;
	_velocityTextureSrc = 0;
	_levelSetTextureSrc = 0;
	_simulationRunning = false;
	_oneStepForward = false;
	_reinitLevelSetNow = false;

	_reinitLevelSetInterval = 10;
	_framesSinceLastReinitCounter = 0;

	_fontRenderer = new FontRenderer();
	_fontRenderer->loadFont("fonts/LucidaSansTypewriter.ttf", "Lucida", 12);
	_fontRenderer->setWindowSize((float)_windowWidth, (float)_windowHeight);
	Text* numberText = _fontRenderer->createText("13.2346234", "Lucida", 200);
	numberText->setPosition(100, 100);
	numberText->setColor(0.f, 0.f, 0.f);
	_numberTextBox = new TextBox(numberText);
	TextBox::setWindowSize((float)_windowWidth, (float)_windowHeight);
	
	glBlendEquation(GL_FUNC_ADD);
	glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);

	/*
	PCGSolverTest* test = new PCGSolverTest();
	test->test3();
	delete test;
	*/

	_inputHandler = inputHandler;
	_inputHandler->addListener((KeyEventListener*) this);
	_inputHandler->addListener((WindowEventListener*) this);
	_inputHandler->addListener((MouseEventListener*) this);

	_clearDataProgram = new GLProgram("Clear Data Program");
	_clearDataProgram->attach(GL_COMPUTE_SHADER, "shaders/cleardata.glsl");
	_clearDataProgram->link();
	_clearDataProgramUL.n = _clearDataProgram->getUniformLocation("n");
	_clearDataProgramUL.value = _clearDataProgram->getUniformLocation("value");

	_clearImageR32iProgram = new GLProgram("Clear Image r32i Program");
	_clearImageR32iProgram->attach(GL_COMPUTE_SHADER, "shaders/clearimager32i.glsl");
	_clearImageR32iProgram->link();
	_clearImageR32iProgramUL.bounds = _clearImageR32iProgram->getUniformLocation("bounds");
	_clearImageR32iProgramUL.value = _clearImageR32iProgram->getUniformLocation("value");

	_clearImageRgba32fProgram = new GLProgram("Clear Image rgba32f Program");
	_clearImageRgba32fProgram->attach(GL_COMPUTE_SHADER, "shaders/clearimagergba32f.glsl");
	_clearImageRgba32fProgram->link();
	_clearImageRgba32fProgramUL.bounds = _clearImageRgba32fProgram->getUniformLocation("bounds");
	_clearImageRgba32fProgramUL.value = _clearImageRgba32fProgram->getUniformLocation("value");

	_pcgSolver = new PCGSolver();

	for (unsigned i=0;i<2;i++)
	{
		_levelSetTexture[i] = GLTexture::newImmutable(GL_TEXTURE_2D, GL_R32F, GRID_SIZE_X+1, GRID_SIZE_Y+1);
		_levelSetTexture[i]->setFilters(GL_LINEAR, GL_LINEAR);
		_levelSetTexture[i]->setWrap(GL_CLAMP_TO_EDGE);
		//float levelSetBorder[] = {1.0f, 0.0f, 0.0f, 0.0f};
		//glTexParameterfv(GL_TEXTURE_2D, GL_TEXTURE_BORDER_COLOR, levelSetBorder);
	}
	//_levelSetFramebuffer = new GLFramebuffer();
	//_levelSetFramebuffer->attach(_levelSetTexture, GL_COLOR_ATTACHMENT0);
	//_levelSetFramebuffer->setDrawBuffer(GL_COLOR_ATTACHMENT0);

	_levelSetStateTexture = GLTexture::newImmutable(GL_TEXTURE_2D, GL_R32UI, GRID_SIZE_X + 1, GRID_SIZE_Y + 1);

	for (unsigned i = 0; i < 2; i++)
		_levelSetActiveNodesBuffer[i] = GLBuffer::newMutable(GL_SHADER_STORAGE_BUFFER, sizeof(GLuint) * 2 * NUM_NODES, 0, GL_STREAM_COPY);

	_levelSetCounterBuffer = GLBuffer::newMutable(GL_ATOMIC_COUNTER_BUFFER, sizeof(GLuint) * 20, 0, GL_STREAM_COPY);
	_levelSetNewValuesBuffer = GLBuffer::newMutable(GL_SHADER_STORAGE_BUFFER, sizeof(GLfloat)*NUM_NODES, 0, GL_STREAM_COPY);

	for (unsigned i=0; i<2; i++)
	{
		_velocityTexture[i] = GLTexture::newImmutable(GL_TEXTURE_2D, GL_RG32F, GRID_SIZE_X+1, GRID_SIZE_Y+1);
		_velocityTexture[i]->setFilters(GL_LINEAR, GL_LINEAR);
		_velocityTexture[i]->setWrap(GL_CLAMP_TO_EDGE);
	}

	for (unsigned i=0;i<2;i++)
	{
		_jumpFloodingTexture[i] = GLTexture::newImmutable(GL_TEXTURE_2D, GL_RGBA32F, GRID_SIZE_X+1, GRID_SIZE_Y+1);
		_jumpFloodingTexture[i]->setFilters(GL_NEAREST, GL_NEAREST);
		_jumpFloodingTexture[i]->setWrap(GL_REPEAT);
	}

	_jumpFloodingFramebuffer = new GLFramebuffer();
	for (unsigned i=0;i<2;i++)
		_jumpFloodingFramebuffer->attach(_jumpFloodingTexture[i], GL_COLOR_ATTACHMENT0+i);

	_nodeIndexTexture = GLTexture::newImmutable(GL_TEXTURE_2D, GL_R32I, GRID_SIZE_X+1, GRID_SIZE_Y+1);

	_femElementPositionsBuffer = GLBuffer::newMutable(GL_SHADER_STORAGE_BUFFER, sizeof(GLuint)*2*NUM_ELEMENTS, NULL, GL_STREAM_COPY);

	_sparseMatrixBuffer = GLBuffer::newMutable(GL_SHADER_STORAGE_BUFFER, sizeof(GLfloat)*9*NUM_NODES, NULL, GL_STREAM_COPY);
	_bVectorBuffer = GLBuffer::newMutable(GL_SHADER_STORAGE_BUFFER, sizeof(GLfloat)*NUM_NODES, NULL, GL_STREAM_COPY);
	_nodePositionBuffer = GLBuffer::newMutable(GL_SHADER_STORAGE_BUFFER, sizeof(GLuint)*2*NUM_NODES, NULL, GL_STREAM_COPY);
	_nodeNeighborhoodBuffer = GLBuffer::newMutable(GL_SHADER_STORAGE_BUFFER, sizeof(GLint)*9*NUM_NODES, NULL, GL_STREAM_COPY);
	_xVectorBuffer = GLBuffer::newMutable(GL_SHADER_STORAGE_BUFFER, sizeof(GLfloat)*NUM_NODES, NULL, GL_STREAM_COPY);

	glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 0, _xVectorBuffer->name);
	_clearDataProgram->use();
	glUniform1ui(_clearDataProgramUL.n, NUM_NODES);
	glUniform1f(_clearDataProgramUL.value, 0.0f);
	glDispatchCompute((unsigned)ceil(NUM_NODES/256.f),1,1);

	_levelSetProgram = new GLProgram("Level Set Program");
	_levelSetProgram->attach(GL_VERTEX_SHADER, "shaders/levelset.vert");
	_levelSetProgram->attach(GL_FRAGMENT_SHADER, "shaders/levelset.frag");
	_levelSetProgram->link();
	_levelSetProgramUL.viewProjMatrix = _levelSetProgram->getUniformLocation("viewProjMatrix");

	_gridProgram = new GLProgram("Grid Program");
	_gridProgram->attach(GL_VERTEX_SHADER, "shaders/grid.vert");
	_gridProgram->attach(GL_FRAGMENT_SHADER, "shaders/grid.frag");
	_gridProgram->link();
	_gridProgramUL.viewProjMatrix = _gridProgram->getUniformLocation("viewProjMatrix");
	_gridProgramUL.color = _gridProgram->getUniformLocation("color");

	_velocityProgram = new GLProgram("Velocity Program");
	_velocityProgram->attach(GL_VERTEX_SHADER, "shaders/velocity.vert");
	_velocityProgram->attach(GL_GEOMETRY_SHADER, "shaders/velocity.geom");
	_velocityProgram->attach(GL_FRAGMENT_SHADER, "shaders/velocity.frag");
	_velocityProgram->link();
	_velocityProgramUL.viewProjMatrix = _velocityProgram->getUniformLocation("viewProjMatrix");
	_velocityProgramUL.color = _velocityProgram->getUniformLocation("color");

	_fluidNodesProgram = new GLProgram("Fluid Nodes Program");
	_fluidNodesProgram->attach(GL_VERTEX_SHADER, "shaders/fluidnodes.vert");
	_fluidNodesProgram->attach(GL_FRAGMENT_SHADER, "shaders/fluidnodes.frag");
	_fluidNodesProgram->link();
	_fluidNodesProgramUL.viewProjMatrix = _fluidNodesProgram->getUniformLocation("viewProjMatrix");
	_fluidNodesProgramUL.color = _fluidNodesProgram->getUniformLocation("color");

	_markFluidNodesProgram = new GLProgram("Mark Fluid Nodes Program");
	_markFluidNodesProgram->attach(GL_COMPUTE_SHADER, "shaders/markfluidnodes.glsl");
	_markFluidNodesProgram->link();

	_initJumpFloodingProgram = new GLProgram("Init Jump Flooding Program");
	_initJumpFloodingProgram->attach(GL_COMPUTE_SHADER, "shaders/initjumpflooding.glsl");
	_initJumpFloodingProgram->link();

	_jumpFloodingProgram = new GLProgram("Jump Flooding Program");
	_jumpFloodingProgram->attach(GL_COMPUTE_SHADER, "shaders/jumpflooding.glsl");
	_jumpFloodingProgram->link();
	_jumpFloodingProgramUL.k = _jumpFloodingProgram->getUniformLocation("k");

	_jumpFloodingFinalProgram = new GLProgram("Jump Flooding Final Program");
	_jumpFloodingFinalProgram->attach(GL_COMPUTE_SHADER, "shaders/jumpfloodingfinal.glsl");
	_jumpFloodingFinalProgram->link();

	_debugTextureProgram = new GLProgram("Debug Texture Program");
	_debugTextureProgram->attach(GL_VERTEX_SHADER, "shaders/debugtexture.vert");
	_debugTextureProgram->attach(GL_FRAGMENT_SHADER, "shaders/debugtexture.frag");
	_debugTextureProgram->link();
	_debugTextureProgramUL.viewProjMatrix = _debugTextureProgram->getUniformLocation("viewProjMatrix");
	
	_advectAndForceProgram = new GLProgram("Advect and Force Program");
	_advectAndForceProgram->attach(GL_COMPUTE_SHADER, "shaders/advectAndForce.glsl");
	_advectAndForceProgram->link();
	_advectAndForceProgramUL.timeStep = _advectAndForceProgram->getUniformLocation("timeStep");

	_solidBoundaryProgram = new GLProgram("Solid Boundary Program");
	_solidBoundaryProgram->attach(GL_COMPUTE_SHADER, "shaders/solidboundary.glsl");
	_solidBoundaryProgram->link();
	
	_enumerateNodesProgram = new GLProgram("Enumerate Nodes Program");
	_enumerateNodesProgram->attach(GL_COMPUTE_SHADER, "shaders/enumeratenodes.glsl");
	_enumerateNodesProgram->link();

	_assembleFemSystemProgram = new GLProgram("Assemble FEM System Program");
	_assembleFemSystemProgram->attach(GL_COMPUTE_SHADER, "shaders/assembleFemSystem.glsl");
	_assembleFemSystemProgram->link();
	_assembleFemSystemProgramUL.numFemElements = _assembleFemSystemProgram->getUniformLocation("numFemElements");

	_generateNodeNeighborhoodProgram = new GLProgram("Assemble FEM System Program");
	_generateNodeNeighborhoodProgram->attach(GL_COMPUTE_SHADER, "shaders/generateNodeNeighborhood.glsl");
	_generateNodeNeighborhoodProgram->link();
	_generateNodeNeighborhoodProgramUL.numFemNodes = _generateNodeNeighborhoodProgram->getUniformLocation("numFemNodes");

	_projectPressureProgram = new GLProgram("Project Pressure Program");
	_projectPressureProgram->attach(GL_COMPUTE_SHADER, "shaders/projectpressure.glsl");
	_projectPressureProgram->link();
	_projectPressureProgramUL.numFemNodes = _projectPressureProgram->getUniformLocation("numFemNodes");

	_levelSetMarkBoundaryProgram = new GLProgram("Level Set Mark Boundary Program");
	_levelSetMarkBoundaryProgram->attach(GL_COMPUTE_SHADER, "shaders/levelsetMarkBoundary.glsl");
	_levelSetMarkBoundaryProgram->link();

	_levelSetMarkUnmarkedNeighborsProgram = new GLProgram("Level Set Mark Unmarked Neighbors Program");
	_levelSetMarkUnmarkedNeighborsProgram->attach(GL_COMPUTE_SHADER, "shaders/levelsetMarkUnmarkedNeighbors.glsl");
	_levelSetMarkUnmarkedNeighborsProgram->link();
	_levelSetMarkUnmarkedNeighborsProgramUL.numFinishedNodes = _levelSetMarkUnmarkedNeighborsProgram->getUniformLocation("numFinishedNodes");

	_levelSetComputeEikonalProgram = new GLProgram("Level Set Compute Eikonal Program");
	_levelSetComputeEikonalProgram->attach(GL_COMPUTE_SHADER, "shaders/levelsetComputeEikonal.glsl");
	_levelSetComputeEikonalProgram->link();
	_levelSetComputeEikonalProgramUL.numActiveNodes = _levelSetComputeEikonalProgram->getUniformLocation("numActiveNodes");
	_levelSetComputeEikonalProgramUL.findNextActiveNodes = _levelSetComputeEikonalProgram->getUniformLocation("findNextActiveNodes");

	_levelSetWriteNewValuesProgram = new GLProgram("Level Set Write New Values Program");
	_levelSetWriteNewValuesProgram->attach(GL_COMPUTE_SHADER, "shaders/levelsetWriteNewValues.glsl");
	_levelSetWriteNewValuesProgram->link();
	_levelSetWriteNewValuesProgramUL.numActiveNodes = _levelSetWriteNewValuesProgram->getUniformLocation("numActiveNodes");

	_levelSetCopyOldValuesProgram = new GLProgram("Level Set Copy Old Values Program");
	_levelSetCopyOldValuesProgram->attach(GL_COMPUTE_SHADER, "shaders/levelsetCopyOldValues.glsl");
	_levelSetCopyOldValuesProgram->link();

	initLevelSetTexture();
	initLevelSetTransferTexture();
	initVelocityTexture();
	initZeroOneQuad();
	initGrid();
	initPointGrid();

	_counterBuffer = GLBuffer::newMutable(GL_ATOMIC_COUNTER_BUFFER, sizeof(GLuint)*20, NULL, GL_STREAM_READ);

	_cameraCenter = glm::vec2(0.5f, 0.5f);
	updateViewMatrix();
	updateProjMatrix();

	_interval = 0.0166666f;

	extrapolateVelocities();
	enforceSolidBoundary();

	return 1;
}

void printFloatBuffer(GLBuffer* buffer, unsigned n)
{
	glMemoryBarrier(GL_ALL_BARRIER_BITS);

	glBindBuffer(GL_SHADER_STORAGE_BUFFER, buffer->name);
	GLfloat* temp = new GLfloat[n];
	glGetBufferSubData(GL_SHADER_STORAGE_BUFFER, 0, sizeof(GLfloat)*n, temp);
	for (unsigned i=0;i<n;i++)
		printf("%f\n",temp[i]);
	printf("\n");
	delete[] temp;
}



void App::buildFEMSystem()
{
	glMemoryBarrier(GL_SHADER_IMAGE_ACCESS_BARRIER_BIT);

	glBindImageTexture(0, _nodeIndexTexture->name, 0, GL_FALSE, 0, GL_WRITE_ONLY, GL_R32I);
	_clearImageR32iProgram->use();
	glUniform2ui(_clearImageR32iProgramUL.bounds, GRID_SIZE_X+1, GRID_SIZE_Y+1);
	glUniform1i(_clearImageR32iProgramUL.value, -1);
	glDispatchCompute((unsigned)ceil((GRID_SIZE_X+1)/16.f), (unsigned)ceil((GRID_SIZE_Y+1)/16.f), 1);

	glMemoryBarrier(GL_ATOMIC_COUNTER_BARRIER_BIT | GL_SHADER_IMAGE_ACCESS_BARRIER_BIT | GL_BUFFER_UPDATE_BARRIER_BIT);

	GLuint zero=0u;
	glBindBufferBase(GL_ATOMIC_COUNTER_BUFFER, 0, _counterBuffer->name);
	glClearBufferData(GL_ATOMIC_COUNTER_BUFFER, GL_R32UI, GL_RED_INTEGER, GL_UNSIGNED_INT, &zero);

	glMemoryBarrier(GL_BUFFER_UPDATE_BARRIER_BIT | GL_TEXTURE_FETCH_BARRIER_BIT);

	glBindImageTexture(0, _levelSetTexture[_levelSetTextureSrc]->name, 0, GL_FALSE, 0, GL_READ_ONLY, GL_R32F);
	glActiveTexture(GL_TEXTURE0);
	glBindTexture(GL_TEXTURE_2D, _levelSetTexture[_levelSetTextureSrc]->name);
	glBindImageTexture(1, _nodeIndexTexture->name, 0, GL_FALSE, 0, GL_READ_WRITE, GL_R32I);

	glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 0, _femElementPositionsBuffer->name);

	// Mark nodes adjacent to fluid
	_markFluidNodesProgram->use();
	glDispatchCompute((unsigned) ceil(GRID_SIZE_X/16.f), (unsigned) ceil(GRID_SIZE_Y/16.f), 1);

	glMemoryBarrier(GL_TEXTURE_FETCH_BARRIER_BIT | GL_ATOMIC_COUNTER_BARRIER_BIT | GL_SHADER_IMAGE_ACCESS_BARRIER_BIT | GL_SHADER_STORAGE_BARRIER_BIT);

	GLuint numFemElements;
	glGetBufferSubData(GL_ATOMIC_COUNTER_BUFFER, 4, sizeof(GLuint), &numFemElements);

	//printf("Num elements: %u\n", numFemElements);

	/*GLint* elemPositions = new GLint[2*numFemElements];
	glGetBufferSubData(GL_SHADER_STORAGE_BUFFER, 0, sizeof(GLint)*2*numFemElements, elemPositions);
	printf("Element positions:\n");
	for (unsigned i=0; i<numFemElements; i++)
		printf("Elem %u: %u, %u\n", i, elemPositions[2*i], elemPositions[2*i+1]);
	printf("\n");
	delete[] elemPositions;*/


	// enumerate nodes that are rows in the FEM system
	glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 3, _nodePositionBuffer->name);
	_enumerateNodesProgram->use();
	glDispatchCompute((unsigned) ceil((GRID_SIZE_X+1)/16.f), (unsigned) ceil((GRID_SIZE_Y+1)/16.f), 1);

	glMemoryBarrier(GL_SHADER_IMAGE_ACCESS_BARRIER_BIT | GL_SHADER_STORAGE_BARRIER_BIT);
	
	GLuint numFemNodes;
	glGetBufferSubData(GL_ATOMIC_COUNTER_BUFFER, 0, sizeof(GLuint), &numFemNodes);

	/*GLuint* nodePositions = new GLuint[2*numFemNodes];
	glGetBufferSubData(GL_SHADER_STORAGE_BUFFER, 0, sizeof(GLuint)*2*numFemNodes, nodePositions);
	for (unsigned i=0; i<numFemNodes; i++)
		printf("Node %u: %u, %u\n", i, nodePositions[2*i], nodePositions[2*i+1]);
	printf("\n");
	delete[] nodePositions;*/

	// generate neighborhood information
	glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 4, _nodeNeighborhoodBuffer->name);
	_generateNodeNeighborhoodProgram->use();
	glUniform1ui(_generateNodeNeighborhoodProgramUL.numFemNodes, numFemNodes);
	glDispatchCompute((unsigned)ceil(numFemNodes/128.f),1,1);

	glMemoryBarrier(GL_SHADER_STORAGE_BARRIER_BIT);

	/*GLint* neighborhood = new GLint[9*numFemNodes];
	glGetBufferSubData(GL_SHADER_STORAGE_BUFFER, 0, sizeof(GLint)*9*numFemNodes, neighborhood);
	for (unsigned i=0; i<numFemNodes; i++)
	{
		printf("Node %u: ", i);
		for (unsigned j=0;j<9;j++)
		{
			int val = neighborhood[9*i+j];
			if (val == -1)
				printf("x ");
			else
				printf("%i ", val);
		}
		printf("\n");
	}
	delete[] neighborhood;*/

	// assemble system matrix and rhs
	float floatZero = 0.f;

	glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 0, _sparseMatrixBuffer->name);
	_clearDataProgram->use();
	glUniform1ui(_clearDataProgramUL.n, 9*numFemNodes);
	glUniform1f(_clearDataProgramUL.value, 0.0f);
	glDispatchCompute((unsigned)ceil(9.f*numFemNodes/256.f),1,1);

	glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 0, _bVectorBuffer->name);
	glUniform1ui(_clearDataProgramUL.n, numFemNodes);
	glDispatchCompute((unsigned)ceil(numFemNodes/256.f),1,1);

	glMemoryBarrier(GL_SHADER_STORAGE_BARRIER_BIT | GL_SHADER_IMAGE_ACCESS_BARRIER_BIT);

	glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 0, _femElementPositionsBuffer->name);
	glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 1, _sparseMatrixBuffer->name);
	glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 2, _bVectorBuffer->name);

	glBindImageTexture(2, _velocityTexture[_velocityTextureSrc]->name, 0, GL_FALSE, 0, GL_READ_ONLY, GL_RG32F);

	_assembleFemSystemProgram->use();
	glUniform1ui(_assembleFemSystemProgramUL.numFemElements, numFemElements);
	glDispatchCompute((unsigned)ceil(numFemElements/128.f), 1, 1);

	glMemoryBarrier(GL_ALL_BARRIER_BITS);

	_numFemNodes = numFemNodes;
	_numFemElements = numFemElements;

	
	//glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 1, 0);
	//glBindBuffer(GL_SHADER_STORAGE_BUFFER, 0);
	
	/*
	glBindBuffer(GL_SHADER_STORAGE_BUFFER, _sparseMatrixBuffer->name);
	GLfloat* sparseMat = new GLfloat[9*numFemNodes];
	glGetBufferSubData(GL_SHADER_STORAGE_BUFFER, 0, sizeof(GLfloat)*9*numFemNodes, sparseMat);
	printf("Sparse Matrix:\n");
	for (unsigned i=0; i<numFemNodes; i++)
	{
		printf("Node %u: ", i);
		for (unsigned j=0;j<9;j++)
		{
			GLfloat val = sparseMat[9*i+j];
			printf("%f ", val);
		}
		printf("\n");
	}
	printf("\n");
	delete[] sparseMat;

	glBindBuffer(GL_SHADER_STORAGE_BUFFER, _bVectorBuffer->name);

	GLfloat* bVector = new GLfloat[numFemNodes];
	glGetBufferSubData(GL_SHADER_STORAGE_BUFFER, 0, sizeof(GLfloat)*numFemNodes, bVector);
	printf("bVector:\n");
	for(unsigned i=0;i<numFemNodes;i++)
		printf("%f ", bVector[i]);
	printf("\n");
	delete[] bVector;
	*/
}

void App::advectAndApplyForces()
{
	glActiveTexture(GL_TEXTURE0);
	glBindTexture(GL_TEXTURE_2D, _velocityTexture[_velocityTextureSrc]->name);
	glActiveTexture(GL_TEXTURE1);
	glBindTexture(GL_TEXTURE_2D, _levelSetTexture[_levelSetTextureSrc]->name);
	glBindImageTexture(0, _velocityTexture[1-_velocityTextureSrc]->name, 0, GL_FALSE, 0, GL_WRITE_ONLY, GL_RG32F);
	glBindImageTexture(1, _levelSetTexture[1-_levelSetTextureSrc]->name, 0, GL_FALSE, 0, GL_WRITE_ONLY, GL_R32F);

	_advectAndForceProgram->use();
	glUniform1f(_advectAndForceProgramUL.timeStep, (float)_interval);

	glMemoryBarrier(GL_TEXTURE_FETCH_BARRIER_BIT | GL_SHADER_IMAGE_ACCESS_BARRIER_BIT);

	glDispatchCompute((unsigned) ceil((GRID_SIZE_X+1)/16.f), (unsigned) ceil((GRID_SIZE_Y+1)/16.f), 1);

	_levelSetTextureSrc = 1-_levelSetTextureSrc;
	_velocityTextureSrc = 1-_velocityTextureSrc;
}

void App::enforceSolidBoundary()
{
	glBindImageTexture(0, _velocityTexture[_velocityTextureSrc]->name, 0, GL_FALSE, 0, GL_READ_WRITE, GL_RG32F);
	_solidBoundaryProgram->use();

	glMemoryBarrier(GL_SHADER_IMAGE_ACCESS_BARRIER_BIT);

	glDispatchCompute((unsigned) ceil((GRID_SIZE_X+1)/16.f), (unsigned) ceil((GRID_SIZE_Y+1)/16.f), 1);
}

void App::extrapolateVelocities()
{
	/*
	GLuint query;
	glGenQueries(1, &query);
	glBeginQuery(GL_TIME_ELAPSED, query);
	*/

	glMemoryBarrier(GL_SHADER_IMAGE_ACCESS_BARRIER_BIT);

	int longestSide = glm::max(GRID_SIZE_X,GRID_SIZE_Y);
	float maxDist = 4.f * longestSide * longestSide;
	float fillValue[4] = {maxDist, maxDist, maxDist, 1.f};

	glBindImageTexture(0, _jumpFloodingTexture[0]->name, 0, GL_FALSE, 0, GL_WRITE_ONLY, GL_RGBA32F);
	_clearImageRgba32fProgram->use();
	glUniform2ui(_clearImageRgba32fProgramUL.bounds, GRID_SIZE_X+1, GRID_SIZE_Y+1);
	glUniform4f(_clearImageRgba32fProgramUL.value, maxDist, maxDist, maxDist, 1.f);
	glDispatchCompute((unsigned)ceil((GRID_SIZE_X+1)/16.f), (unsigned)ceil((GRID_SIZE_Y+1)/16.f), 1);

	glBindImageTexture(0, _levelSetTexture[_levelSetTextureSrc]->name, 0, GL_FALSE, 0, GL_READ_ONLY, GL_R32F);
	glBindImageTexture(1, _jumpFloodingTexture[0]->name, 0, GL_FALSE, 0, GL_WRITE_ONLY, GL_RGBA32F);

	//_markFluidNodesProgram->use();
	//glDispatchCompute((unsigned) ceil(GRID_SIZE_X/16.f), (unsigned) ceil(GRID_SIZE_Y/16.f), 1);

	glMemoryBarrier(GL_SHADER_IMAGE_ACCESS_BARRIER_BIT);
	
	_initJumpFloodingProgram->use();
	glDispatchCompute((unsigned) ceil((GRID_SIZE_X+1)/16.f), (unsigned) ceil((GRID_SIZE_Y+1)/16.f), 1);

	glMemoryBarrier(GL_TEXTURE_FETCH_BARRIER_BIT | GL_SHADER_IMAGE_ACCESS_BARRIER_BIT);

	_jumpFloodingSrc = 0;

	int src=0;
	glm::uvec2 k(GRID_SIZE_X+1, GRID_SIZE_Y+1);
	k = glm::uvec2((unsigned) ceil(k.x/2.f), (unsigned) ceil(k.y/2.f));

	_jumpFloodingProgram->use();
	glActiveTexture(GL_TEXTURE0);
	do {
		int dst = 1-src;
		glBindTexture(GL_TEXTURE_2D, _jumpFloodingTexture[src]->name);
		glBindImageTexture(0, _jumpFloodingTexture[dst]->name, 0, GL_FALSE, 0, GL_WRITE_ONLY, GL_RGBA32F);
		glUniform2ui(_jumpFloodingProgramUL.k, k.x, k.y);

		glMemoryBarrier(GL_TEXTURE_FETCH_BARRIER_BIT | GL_SHADER_IMAGE_ACCESS_BARRIER_BIT);

		glDispatchCompute((unsigned) ceil((GRID_SIZE_X+1)/16.f), (unsigned) ceil((GRID_SIZE_Y+1)/16.f), 1);

		src = 1-src;
		k = glm::uvec2((unsigned) ceil(k.x/2.f), (unsigned) ceil(k.y/2.f));
	} while(k.x > 1 || k.y > 1);

	_jumpFloodingSrc = src;

	// _jumpFloodingTexture[src] contains the final distSq values and coordinates of the closest fluid point
	glBindImageTexture(0, _jumpFloodingTexture[src]->name, 0, GL_FALSE, 0, GL_READ_ONLY, GL_RGBA32F);
	glBindImageTexture(1, _velocityTexture[_velocityTextureSrc]->name, 0, GL_FALSE, 0, GL_READ_WRITE, GL_RG32F);

	glMemoryBarrier(GL_SHADER_IMAGE_ACCESS_BARRIER_BIT);

	_jumpFloodingFinalProgram->use();
	glDispatchCompute((unsigned) ceil((GRID_SIZE_X+1)/16.f), (unsigned) ceil((GRID_SIZE_Y+1)/16.f), 1);

	/*
	glEndQuery(GL_TIME_ELAPSED);
	GLuint result;
	glGetQueryObjectuiv(query, GL_QUERY_RESULT, &result);
	glDeleteQueries(1, &query);
	printf("Velocity Extrapolation: %f ms\n", result / 1000000.f);
	*/
}

void App::initVelocityTexture()
{
	glm::vec2* velocityData = new glm::vec2[NUM_NODES];

	for (unsigned y=0; y<=GRID_SIZE_Y; y++)
	{
		for (unsigned x=0; x<=GRID_SIZE_X; x++)
		{
			// Radial outwards velocity pattern
			//velocityData[NODE_INDEX(x,y)] = glm::vec2((1.3f/GRID_SIZE_X)*(x-GRID_SIZE_X/2.f), (1.3f/GRID_SIZE_Y)*(y-GRID_SIZE_Y/2.f));

			// Vertical stripe of movement
			//velocityData[NODE_INDEX(x,y)] = glm::vec2(1.f*exp(-pow((x-0.5f*GRID_SIZE_X)/(0.05f*GRID_SIZE_X),2)), 0.f);
			
			// Zero everything
			velocityData[NODE_INDEX(x,y)] = glm::vec2(0.f,0.f);

		}
	}

	for (unsigned x=0; x < GRID_SIZE_X/2; x++)
	{
		velocityData[NODE_INDEX(x,(unsigned)(GRID_SIZE_Y*2.f/3.f))] = glm::vec2(10.f, 0.f);
		velocityData[NODE_INDEX(x,(unsigned)(GRID_SIZE_Y*2.f/3.f)+1)] = glm::vec2(10.f, 0.f);
	}

	for (unsigned x=GRID_SIZE_X/2; x<=GRID_SIZE_X; x++)
	{
		velocityData[NODE_INDEX(x,(unsigned)(GRID_SIZE_Y*2.f/3.f))] = glm::vec2(-10.f, 0.f);
		velocityData[NODE_INDEX(x,(unsigned)(GRID_SIZE_Y*2.f/3.f)+1)] = glm::vec2(-10.f, 0.f);
	}

	_velocityTexture[_velocityTextureSrc]->setContents(GL_RG, GL_FLOAT, &velocityData[0][0]);

	delete[] velocityData;

	
}

void App::initZeroOneQuad()
{
	// [0,0]x[1,1] quad
	float positionData[16] =
	{
		0.f,0.f,0.f,1.f,
		1.f,0.f,0.f,1.f,
		0.f,1.f,0.f,1.f,
		1.f,1.f,0.f,1.f
	};
	_zeroOnePositionBuffer = GLBuffer::newMutable(GL_ARRAY_BUFFER, sizeof(float)*16, positionData, GL_STATIC_DRAW);
	_zeroOneVertexArray = new GLVertexArray();
	_zeroOneVertexArray->attach(_zeroOnePositionBuffer, 0, 4, GL_FLOAT);
}

void App::initPointGrid()
{
	// one point per node in [0,0]x[1,1]
	unsigned size = 4*NUM_NODES;
	float* positionData = new float[size];

	for (unsigned y=0; y<=GRID_SIZE_Y; y++)
	{
		for (unsigned x=0; x<=GRID_SIZE_X; x++)
		{
			positionData[NODE_INDEX(x,y)*4+0] = (float)x / (GRID_SIZE_X);
			positionData[NODE_INDEX(x,y)*4+1] = (float)y / (GRID_SIZE_Y);
			positionData[NODE_INDEX(x,y)*4+2] = 0.f;
			positionData[NODE_INDEX(x,y)*4+3] = 1.f;
		}
	}

	_pointGridBuffer = GLBuffer::newMutable(GL_ARRAY_BUFFER, sizeof(float)*size, positionData, GL_STATIC_DRAW);
	_pointGridVertexArray = new GLVertexArray();
	_pointGridVertexArray->attach(_pointGridBuffer, 0, 4, GL_FLOAT);

	delete[] positionData;
}

void App::initGrid()
{
	// line grid to visualize elements
	unsigned size = 4*2*(GRID_SIZE_X+1+GRID_SIZE_Y+1);
	float* positionData = new float[size];

	for (unsigned i=0;i<GRID_SIZE_X+1;i++)
	{
		positionData[8*i+0] = (float)i/GRID_SIZE_X;
		positionData[8*i+1] = 0.f;
		positionData[8*i+2] = 0.f;
		positionData[8*i+3] = 1.f;
		positionData[8*i+4] = (float)i/GRID_SIZE_X;
		positionData[8*i+5] = 1.f;
		positionData[8*i+6] = 0.f;
		positionData[8*i+7] = 1.f;
	}

	unsigned offset = 8*(GRID_SIZE_X+1);
	for (unsigned i=0;i<GRID_SIZE_Y+1;i++)
	{
		positionData[offset+8*i+0] = 0.f;
		positionData[offset+8*i+1] = (float)i/GRID_SIZE_Y;
		positionData[offset+8*i+2] = 0.f;
		positionData[offset+8*i+3] = 1.f;
		positionData[offset+8*i+4] = 1.f;
		positionData[offset+8*i+5] = (float)i/GRID_SIZE_Y;
		positionData[offset+8*i+6] = 0.f;
		positionData[offset+8*i+7] = 1.f;
	}

	_gridBuffer = GLBuffer::newMutable(GL_ARRAY_BUFFER, sizeof(float)*size, positionData, GL_STATIC_DRAW);
	_gridVertexArray = new GLVertexArray();
	_gridVertexArray->attach(_gridBuffer, 0, 4, GL_FLOAT);

	delete[] positionData;
}

void App::initLevelSetTexture()
{
	// disk of water in upper half
	float* levelSet = new float[NUM_NODES];

	float center[2] = {GRID_SIZE_X, GRID_SIZE_Y * (1.f/6.f)};
	float radius = GRID_SIZE_Y * (1.f/3.f);

	glm::vec2 minPoint(GRID_SIZE_X*0.5f-0.5f, GRID_SIZE_Y*0.5f-0.5f);
	glm::vec2 maxPoint(GRID_SIZE_X*0.5f+1.5f, GRID_SIZE_Y*0.5f+1.5f);

	for (unsigned y=0; y<=GRID_SIZE_Y; y++)
	{
		for (unsigned x=0; x<=GRID_SIZE_X; x++)
		{
			float level = sqrt((x-center[0]) * (x-center[0]) + (y-center[1])*(y-center[1])) - radius;
			levelSet[NODE_INDEX(x,y)] = level;
		}
	}
	
	/*
	levelSet[NODE_INDEX(GRID_SIZE_X/2, GRID_SIZE_Y/2)] = -0.5f;
	levelSet[NODE_INDEX(GRID_SIZE_X/2+1, GRID_SIZE_Y/2)] = -0.5f;
	levelSet[NODE_INDEX(GRID_SIZE_X/2+1, GRID_SIZE_Y/2+1)] = -0.5f;
	levelSet[NODE_INDEX(GRID_SIZE_X/2, GRID_SIZE_Y/2+1)] = -0.5f;

	for (unsigned x=0; x < GRID_SIZE_X/2; x++)
	{
		levelSet[NODE_INDEX(x,GRID_SIZE_Y/2)] = minPoint.x-x;
		levelSet[NODE_INDEX(x,GRID_SIZE_Y/2+1)] = minPoint.x-x;
	}

	for (unsigned x = GRID_SIZE_X/2+2; x <= GRID_SIZE_X; x++)
	{
		levelSet[NODE_INDEX(x,GRID_SIZE_Y/2)] = x-maxPoint.x;
		levelSet[NODE_INDEX(x,GRID_SIZE_Y/2+1)] = x-maxPoint.x;
	}

	for (unsigned y=0; y < GRID_SIZE_Y/2; y++)
	{
		levelSet[NODE_INDEX(GRID_SIZE_X/2,y)] = minPoint.y-y;
		levelSet[NODE_INDEX(GRID_SIZE_X/2+1,y)] = minPoint.y-y;
	}

	for (unsigned y=GRID_SIZE_Y/2+2; y <= GRID_SIZE_Y; y++)
	{
		levelSet[NODE_INDEX(GRID_SIZE_X/2,y)] = y-maxPoint.y;
		levelSet[NODE_INDEX(GRID_SIZE_X/2+1,y)] = y-maxPoint.y;
	}

	for (unsigned x=0; x < GRID_SIZE_X/2; x++)
	{
		for (unsigned y=0; y < GRID_SIZE_Y/2; y++)
		{
			levelSet[NODE_INDEX(x,y)] = glm::distance(glm::vec2(x,y), minPoint);
		}
	}

	for (unsigned x=GRID_SIZE_X/2+2; x <= GRID_SIZE_X; x++)
	{
		for (unsigned y=0; y < GRID_SIZE_Y/2; y++)
		{
			levelSet[NODE_INDEX(x,y)] = glm::distance(glm::vec2(x,y), glm::vec2(maxPoint.x, minPoint.y));
		}
	}

	for (unsigned x=0; x < GRID_SIZE_X/2; x++)
	{
		for (unsigned y=GRID_SIZE_Y/2+2; y <= GRID_SIZE_Y; y++)
		{
			levelSet[NODE_INDEX(x,y)] = glm::distance(glm::vec2(x,y), glm::vec2(minPoint.x, maxPoint.y));
		}
	}

	for (unsigned x=GRID_SIZE_X/2+2; x <= GRID_SIZE_X; x++)
	{
		for (unsigned y=GRID_SIZE_Y/2+2; y <= GRID_SIZE_Y; y++)
		{
			levelSet[NODE_INDEX(x,y)] = glm::distance(glm::vec2(x,y), maxPoint);
		}
	}
	*/
	

	_levelSetTexture[_levelSetTextureSrc]->setContents(GL_RED, GL_FLOAT, levelSet);

	delete[] levelSet;
}

void App::initLevelSetTransferTexture()
{
	// blue in water, red outside
	glm::vec3* colors = new glm::vec3[MAX_TRANSFER_LEVEL-MIN_TRANSFER_LEVEL];
	glm::vec3 lightBlue(.8f, .8f, 1.f);
	glm::vec3 darkBlue(.4f, .4f, 1.f);
	glm::vec3 lightRed(1.f, .8f, .8f);
	glm::vec3 darkRed(1.f, .4f, .4f);

	for (unsigned i=0;i<-MIN_TRANSFER_LEVEL;i++)
	{
		float f = i / (-MIN_TRANSFER_LEVEL-1.f);
		colors[i] = lightBlue*(1.f-f) + darkBlue*f;
	}

	for (unsigned i=0;i<MAX_TRANSFER_LEVEL;i++)
	{
		float f = i / (MAX_TRANSFER_LEVEL-1.f);
		colors[i-MIN_TRANSFER_LEVEL] = darkRed*(1.f-f) + lightRed*f;
	}

	_levelSetTransferTexture = GLTexture::newImmutable(GL_TEXTURE_1D, GL_RGB8, MAX_TRANSFER_LEVEL-MIN_TRANSFER_LEVEL, 1, GL_RGB, GL_FLOAT, &(colors[0][0]));
	_levelSetTransferTexture->setFilters(GL_NEAREST, GL_NEAREST);
	_levelSetTransferTexture->setWrap(GL_CLAMP_TO_EDGE);

	delete[] colors;
}

void App::projectPressure()
{
	GLuint buffers[] = {
		_nodeNeighborhoodBuffer->name,
		_xVectorBuffer->name,
		_nodePositionBuffer->name
	};

	glMemoryBarrier(GL_SHADER_STORAGE_BARRIER_BIT | GL_SHADER_IMAGE_ACCESS_BARRIER_BIT);

	glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 0, _nodeNeighborhoodBuffer->name);
	glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 1, _xVectorBuffer->name);
	glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 2, _nodePositionBuffer->name);

	glBindImageTexture(0, _velocityTexture[_velocityTextureSrc]->name, 0, GL_FALSE, 0, GL_READ_WRITE, GL_RG32F);

	_projectPressureProgram->use();
	glUniform1ui(_projectPressureProgramUL.numFemNodes, _numFemNodes);

	glDispatchCompute((unsigned) ceil(_numFemNodes / 512.f), 1, 1);

	glMemoryBarrier(GL_SHADER_STORAGE_BARRIER_BIT | GL_TEXTURE_FETCH_BARRIER_BIT | GL_SHADER_IMAGE_ACCESS_BARRIER_BIT);
}

void App::solveFEMSystem()
{
	

	/*
	GLuint query;
	glGenQueries(1, &query);
	glBeginQuery(GL_TIME_ELAPSED, query);
	*/

	_pcgSolver->set(_sparseMatrixBuffer, _nodeNeighborhoodBuffer, _bVectorBuffer, _numFemNodes);

	glMemoryBarrier(GL_SHADER_STORAGE_BARRIER_BIT);

	glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 0, _xVectorBuffer->name);
	_clearDataProgram->use();
	glUniform1ui(_clearDataProgramUL.n, _numFemNodes);
	glUniform1f(_clearDataProgramUL.value, 0.0f);
	glDispatchCompute((unsigned)ceil(_numFemNodes/256.f),1,1);

	_pcgSolver->solve(_xVectorBuffer, 30);

	/*printf("pressure:\n");
	GLfloat* pressure = new GLfloat[_numFemNodes];
	glBindBuffer(GL_SHADER_STORAGE_BUFFER, _xVectorBuffer->name);
	glGetBufferSubData(GL_SHADER_STORAGE_BUFFER, 0, sizeof(GLfloat)*_numFemNodes, pressure);
	for (unsigned i=0;i<_numFemNodes;i++)
		printf("%f ",pressure[i]);
	printf("\n");*/
	


	/*
	glEndQuery(GL_TIME_ELAPSED);
	GLuint result;
	glGetQueryObjectuiv(query, GL_QUERY_RESULT, &result);
	glDeleteQueries(1, &query);
	printf("Solve FEM System: %f ms\n", result / 1000000.f);
	*/
}

void App::reinitLevelSet()
{
	float highFloat = 10000.0f;
	glClearTexImage(_levelSetTexture[1 - _levelSetTextureSrc]->name, 0, GL_RED, GL_FLOAT, &highFloat);

	glMemoryBarrier(GL_TEXTURE_UPDATE_BARRIER_BIT | GL_BUFFER_UPDATE_BARRIER_BIT);

	unsigned zeroValue = 0u;
	glClearTexImage(_levelSetStateTexture->name, 0, GL_RED_INTEGER, GL_UNSIGNED_INT, &zeroValue);

	glBindBufferBase(GL_ATOMIC_COUNTER_BUFFER, 0, _levelSetCounterBuffer->name);
	glClearBufferData(GL_ATOMIC_COUNTER_BUFFER, GL_R32UI, GL_RED_INTEGER, GL_UNSIGNED_INT, &zeroValue);
	
	glMemoryBarrier(GL_TEXTURE_UPDATE_BARRIER_BIT | GL_BUFFER_UPDATE_BARRIER_BIT | GL_TEXTURE_FETCH_BARRIER_BIT
		| GL_SHADER_IMAGE_ACCESS_BARRIER_BIT | GL_ATOMIC_COUNTER_BARRIER_BIT | GL_SHADER_STORAGE_BARRIER_BIT);

	glActiveTexture(GL_TEXTURE0);
	glBindTexture(GL_TEXTURE_2D, _levelSetTexture[_levelSetTextureSrc]->name);
	glBindImageTexture(0, _levelSetStateTexture->name, 0, GL_FALSE, 0, GL_READ_WRITE, GL_R32UI);
	glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 0, _levelSetActiveNodesBuffer[0]->name);

	_levelSetMarkBoundaryProgram->use();
	glDispatchCompute((unsigned)ceil(GRID_SIZE_X / 16.f), (unsigned)ceil(GRID_SIZE_Y / 16.f), 1);

	glMemoryBarrier(GL_BUFFER_UPDATE_BARRIER_BIT);

	GLuint numFinishedNodes;
	glGetBufferSubData(GL_ATOMIC_COUNTER_BUFFER, 0, sizeof(GLuint), &numFinishedNodes);
	glBufferSubData(GL_ATOMIC_COUNTER_BUFFER, 0, sizeof(GLuint), &zeroValue);

	glMemoryBarrier(GL_BUFFER_UPDATE_BARRIER_BIT | GL_ATOMIC_COUNTER_BARRIER_BIT);

	glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 1, _levelSetActiveNodesBuffer[1]->name);
	glBindImageTexture(1, _levelSetTexture[_levelSetTextureSrc]->name, 0, GL_FALSE, 0, GL_READ_ONLY, GL_R32F); // old level set
	glBindImageTexture(2, _levelSetTexture[1-_levelSetTextureSrc]->name, 0, GL_FALSE, 0, GL_WRITE_ONLY, GL_R32F); // new level set, built from boundary

	_levelSetMarkUnmarkedNeighborsProgram->use();
	glUniform1ui(_levelSetMarkUnmarkedNeighborsProgramUL.numFinishedNodes, numFinishedNodes);

	glDispatchCompute((unsigned)ceil(numFinishedNodes / 64.f), 1, 1);

	glMemoryBarrier(GL_BUFFER_UPDATE_BARRIER_BIT | GL_SHADER_IMAGE_ACCESS_BARRIER_BIT);

	GLuint numActiveNodes;
	glGetBufferSubData(GL_ATOMIC_COUNTER_BUFFER, 0, sizeof(GLuint), &numActiveNodes);
	glBufferSubData(GL_ATOMIC_COUNTER_BUFFER, 0, sizeof(GLuint), &zeroValue);

	unsigned srcIndex = 1;
	unsigned dstIndex = 0;

	_levelSetTextureSrc = 1 - _levelSetTextureSrc;

	glActiveTexture(GL_TEXTURE0);
	glBindTexture(GL_TEXTURE_2D, _levelSetTexture[_levelSetTextureSrc]->name); // new level set, built from boundary

	glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 2, _levelSetNewValuesBuffer->name);

	unsigned LOOPSIZE = 5;
	for (unsigned i = 0; i < LOOPSIZE; i++)
	{
		glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 0, _levelSetActiveNodesBuffer[srcIndex]->name);
		glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 1, _levelSetActiveNodesBuffer[dstIndex]->name);

		_levelSetComputeEikonalProgram->use();
		glUniform1ui(_levelSetComputeEikonalProgramUL.numActiveNodes, numActiveNodes);
		if (i == 0)
		{
			glUniform1ui(_levelSetComputeEikonalProgramUL.findNextActiveNodes, 1);
		}
		else if (i == LOOPSIZE-1)
		{
			glUniform1ui(_levelSetComputeEikonalProgramUL.findNextActiveNodes, 0);
		}

		glMemoryBarrier(GL_BUFFER_UPDATE_BARRIER_BIT | GL_SHADER_IMAGE_ACCESS_BARRIER_BIT | GL_TEXTURE_FETCH_BARRIER_BIT | GL_SHADER_STORAGE_BARRIER_BIT | GL_ATOMIC_COUNTER_BARRIER_BIT);

		glDispatchCompute((unsigned)ceil(numActiveNodes / 64.f), 1, 1);

		_levelSetWriteNewValuesProgram->use();
		glUniform1ui(_levelSetWriteNewValuesProgramUL.numActiveNodes, numActiveNodes);

		glMemoryBarrier(GL_SHADER_IMAGE_ACCESS_BARRIER_BIT | GL_SHADER_STORAGE_BARRIER_BIT);

		glDispatchCompute((unsigned)ceil(numActiveNodes / 64.f), 1, 1);

		if (i < LOOPSIZE - 1)
		{
			glMemoryBarrier(GL_BUFFER_UPDATE_BARRIER_BIT);

			srcIndex = 1 - srcIndex;
			dstIndex = 1 - srcIndex;

			glGetBufferSubData(GL_ATOMIC_COUNTER_BUFFER, 0, sizeof(GLuint), &numActiveNodes);
			glBufferSubData(GL_ATOMIC_COUNTER_BUFFER, 0, sizeof(GLuint), &zeroValue);
		}
		
	}

	glMemoryBarrier(GL_SHADER_IMAGE_ACCESS_BARRIER_BIT);

	_levelSetCopyOldValuesProgram->use();
	glDispatchCompute((unsigned)ceil((GRID_SIZE_X + 1) / 16.f), (unsigned)ceil((GRID_SIZE_Y + 1) / 16.f), 1);

	glMemoryBarrier(GL_ALL_BARRIER_BITS);
	
}

int App::loop()
{
	double time = glfwGetTime();
	_interval = time - _startTime;
	_startTime = time;

	_interval = glm::min(.1, _interval);

	//printf("FPS: %f\n",1.0/_interval);

	if (_windowResized)
	{
		updateProjMatrix();
		_fontRenderer->setWindowSize((float)_windowWidth, (float)_windowHeight);
		TextBox::setWindowSize((float)_windowWidth, (float)_windowHeight);
	}

	if (_cameraMoving)
	{
		double xpos, ypos;
		glfwGetCursorPos(_window, &xpos, &ypos);
		glm::vec2 mousePos = cursorToTex(glm::vec2(xpos, ypos));
		glm::vec2 movement = 2.f * glm::inverse(glm::mat2(_projMatrix)) * (mousePos - _mousePositionStart);
		_cameraCenter = _cameraCenterStart - movement;
		updateViewMatrix();
	}

	if (_reinitLevelSetNow || _reinitLevelSetInterval == _framesSinceLastReinitCounter)
	{
		reinitLevelSet();
		_framesSinceLastReinitCounter = 0;
	}
	_framesSinceLastReinitCounter++;
	_reinitLevelSetNow = false;

	if (_simulationRunning || _oneStepForward)
	{
		if (_oneStepForward)
			_interval = 1.f / 60.f;

		extrapolateVelocities();
		advectAndApplyForces();
		enforceSolidBoundary();
		buildFEMSystem();
		solveFEMSystem();

		projectPressure();
		
	}

	glBindFramebuffer(GL_FRAMEBUFFER, 0);
	glViewport(0, 0, _windowWidth, _windowHeight);

	glClearColor(1.f, 1.f, 1.f, 1.f);
	glClear(GL_COLOR_BUFFER_BIT);

	drawLevelSets();

	if (_drawGrid)
		drawGrid();

	if (_drawVelocities)
		drawVelocities();

	if (_drawFluidNodes)
		drawFluidNodes();

	if (_drawDebugTexture)
		drawDebugTexture();

	if (glfwGetKey(_window, GLFW_KEY_P) == GLFW_PRESS)
	{
		_showOption = SHOW_PRESSURE;
	}
	else if (glfwGetKey(_window, GLFW_KEY_N) == GLFW_PRESS)
	{
		_showOption = SHOW_NEIGHBORHOOD;
	}

	showNumbers();
	_showOption = SHOW_NOTHING;

	_windowResized = false;
	_oneStepForward = false;

	_endTime = glfwGetTime();

	return 1;
}

void App::showNumbers()
{
	if (_showOption == SHOW_NOTHING)
		return;

	double xpos, ypos;
	double xmpos, ympos;
	glfwGetCursorPos(_window, &xmpos, &ympos);
	xpos = (xmpos / _windowWidth)*2.0 - 1.0;
	ypos = 1.0 - 2.0*(ympos / _windowHeight);
	glm::vec4 mousevec(xpos, ypos, 0.0, 1.0);
	glm::vec4 worldpos = glm::inverse(_viewProjMatrix) * mousevec;

	unsigned gx = (unsigned)glm::clamp((int)round(worldpos.x * GRID_SIZE_X), 0, GRID_SIZE_X);
	unsigned gy = (unsigned)glm::clamp((int)round(worldpos.y * GRID_SIZE_Y), 0, GRID_SIZE_Y);

	GLint index;
	glGetTextureSubImage(_nodeIndexTexture->name, 0, gx,gy,0,1,1,1,GL_RED_INTEGER,GL_INT, sizeof(GLint), &index);

	if (index >= 0)
	{
		std::ostringstream ss;
		if (_showOption == SHOW_PRESSURE)
		{
			GLfloat pressure;
			glGetNamedBufferSubData(_xVectorBuffer->name, sizeof(GLfloat)*index, sizeof(GLfloat), &pressure);

			ss << "p: " << pressure;
		}
		else if (_showOption == SHOW_NEIGHBORHOOD)
		{
			GLint neighbors[9];
			glGetNamedBufferSubData(_nodeNeighborhoodBuffer->name, sizeof(GLint) * 9 * index, sizeof(GLint) * 9, neighbors);
			ss << "n: ";
			for (unsigned i = 0; i < 9; i++)
			{
				ss << neighbors[i] << " ";
			}
		}

		glm::vec4 projgridpos((float)gx/GRID_SIZE_X, (float)gy/GRID_SIZE_Y, 0.0f, 1.f);
		projgridpos = _viewProjMatrix * projgridpos;
		
		projgridpos.x = (projgridpos.x + 1.0f)*0.5f*_windowWidth;
		projgridpos.y = (1.0f - projgridpos.y)*0.5f*_windowHeight;

		_numberTextBox->getText()->setText(ss.str());
		_numberTextBox->getText()->setPosition((int)projgridpos.x+5, (int)projgridpos.y-4);

		_numberTextBox->render();
	}
}

void App::updateViewMatrix()
{
	_viewMatrix = glm::mat4(
		1.f, 0.f, 0.f, 0.f,
		0.f, 1.f, 0.f, 0.f,
		0.f, 0.f, 1.f, 0.f,
		-_cameraCenter.x, -_cameraCenter.y, 0.f, 1.f
	);
	_viewProjMatrix = _projMatrix * _viewMatrix;
}

void App::updateProjMatrix()
{
	if (_windowWidth > _windowHeight)
	{
		_projMatrix = glm::mat4(
			_zoomFactor/_windowRatio, 0.f, 0.f, 0.f,
			0.f, _zoomFactor, 0.f, 0.f,
			0.f, 0.f, 1.f, 0.f,
			0.f, 0.f, 0.f, 1.f
		);
	}
	else
	{
		_projMatrix = glm::mat4(
			_zoomFactor, 0.f, 0.f, 0.f,
			0.f, _zoomFactor*_windowRatio, 0.f, 0.f,
			0.f, 0.f, 1.f, 0.f,
			0.f, 0.f, 0.f, 1.f
		);
	}

	_viewProjMatrix = _projMatrix * _viewMatrix;
}

void App::drawGrid()
{
	_gridProgram->use();
	glUniformMatrix4fv(_gridProgramUL.viewProjMatrix, 1, GL_FALSE, &(_viewProjMatrix[0][0]));
	glUniform4f(_gridProgramUL.color, 0.9f, 0.9f, 0.9f, 1.f);

	glBindVertexArray(_gridVertexArray->name);
	glDrawArrays(GL_LINES, 0, 2*(GRID_SIZE_X+1+GRID_SIZE_Y+1));
}

void App::drawLevelSets()
{

	glActiveTexture(GL_TEXTURE0);
	glBindTexture(GL_TEXTURE_2D, _levelSetTexture[_levelSetTextureSrc]->name);
	glActiveTexture(GL_TEXTURE1);
	glBindTexture(GL_TEXTURE_1D, _levelSetTransferTexture->name);

	_levelSetProgram->use();
	glUniformMatrix4fv(_levelSetProgramUL.viewProjMatrix, 1, GL_FALSE, &(_viewProjMatrix[0][0]));

	glBindVertexArray(_zeroOneVertexArray->name);
	glDrawArrays(GL_TRIANGLE_STRIP, 0, 4);
}

void App::drawVelocities()
{
	glActiveTexture(GL_TEXTURE0);
	glBindTexture(GL_TEXTURE_2D, _velocityTexture[_velocityTextureSrc]->name);

	_velocityProgram->use();
	glUniformMatrix4fv(_velocityProgramUL.viewProjMatrix, 1, GL_FALSE, &(_viewProjMatrix[0][0]));
	glUniform4f(_velocityProgramUL.color, 0.0f, 0.0f, 0.0f, 1.f);

	glBindVertexArray(_pointGridVertexArray->name);
	glDrawArrays(GL_POINTS, 0, NUM_NODES);
}

void App::drawFluidNodes()
{
	glActiveTexture(GL_TEXTURE0);
	glBindTexture(GL_TEXTURE_2D, _nodeIndexTexture->name);

	_fluidNodesProgram->use();

	glUniformMatrix4fv(_fluidNodesProgramUL.viewProjMatrix, 1, GL_FALSE, &(_viewProjMatrix[0][0]));
	glUniform4f(_fluidNodesProgramUL.color, 0.0f, 0.5f, 0.0f, 1.f);

	glPointSize(3.f);
	glBindVertexArray(_pointGridVertexArray->name);
	glDrawArrays(GL_POINTS, 0, NUM_NODES);
}

void App::drawDebugTexture()
{
	glActiveTexture(GL_TEXTURE0);
	glBindTexture(GL_TEXTURE_2D, _jumpFloodingTexture[_jumpFloodingSrc]->name);

	_debugTextureProgram->use();

	glUniformMatrix4fv(_debugTextureProgramUL.viewProjMatrix, 1, GL_FALSE, &(_viewProjMatrix[0][0]));

	glBindVertexArray(_zeroOneVertexArray->name);
	glDrawArrays(GL_TRIANGLE_STRIP, 0, 4);
}